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ABSTRACT 

The recently discovered circumbinary planets (Kepler-16 b, Kepler 34-b, Kepler 35-b) represent 
the first direct evidence of the viability of planet formation in circumbinary orbits. We report on 
the results of A-body simulations investigating planetesimal accretion in the Kepler-16 b system, 
focusing on the range of impact velocities under the influence of both stars' gravitational perturbation 
and friction from a putative protoplanetary disk. Our results show that planet formation might be 
effectively inhibited for a large range in semi- major axis (1.75 < ap < 4 AU), suggesting that the 
planetary core must have either migrated from outside 4 AU, or formed in situ very close to its current 
location. 

Subject headings: Planets and satellites: formation, Planets and satellites: dynamical evolution and 
stability 



1. INTRODUCTION 

The discovery of extrasolar planets around main- 
sequence stars is one of the major observational break- 
throughs of the last decade. The size of the planetary 
census, propelled by radial velocity (RV) surveys and 
dedicated missions such as Kepler, has grown to include 
planetary systems where a variety of interesting dynam- 
ical interactions can be observed. Such systems include 
61 exoplanets discovered around stellar binariefl (includ- 
ing both planets orbiting one of the stellar companions 
and circumbinary planets). While for the majority of 
these planets the binarity of the system represents only 
a weak perturbation on the gravitational pull of the cen- 
tral star, a few single-planet systems have been detected 
in binaries with abin S 30 AU (such as HD 41004, Gliese 
86, HD196885 and 7 Cephei), with each planet in a 
circumstellar ("S-type") orbit. Only one multiple sys- 
tem with q bin 5, 100 AU has been found (HD 177830, 
iMeschiari et al.ll2011l) . 

The existence of these systems represents a major 
challenge to the current paradigm of planet formation. 
In fact, a number of simulations attempting to model 
the dynamics of the growth of planetary embryos from 
fcm-sized planetesimals in presence of a binary com- 
panion have hit significant difficulties (among others , 
Marzari fc Schollll2000l iThebault et al1l20ql l200i [20061: 
ThebaultJ 120111: iPaardekooper et al.l 120081: iFragner et al.l 
20111 ). The most important parameter controlling plan- 
etesimal accretion is the mutual encounter velocity; in- 
deed, runaway growth requires it to be less than the 
escape velocity for efficient accretion. The presence of 
the companion can stir up the relative velocity between 
planetesimals, interfering with runaway growth. Rela- 
tive velocity is often excited beyond a fiducial thresh- 
old velocity at which all encounters are erosive, poten- 
tially slowing down planet formation or halting it alto- 
gether. Simulations taking into account a static back- 
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ground gas disk (representing an unperturbed protoplan- 
etary disk at some point in time) initially posited that 
disk-planetesimals interaction induces a phasing of the 
orbits, making the envi ronment more accretion-friendly 
(|Marzari fc Scholl|[2000]) . Nevertheless, if the protoplan- 
ets interact with the gas disk through aerodynamic drag 
alone, the phasing induced by the gas disk is clearly size- 
dependent, and protoplanets with different sizes will col- 
lide with large encounter speeds over the majority of 
the r ange in semi-major axis sampled ()Thebault et al.l 
l2004f ). Finally, a misalignment between the orbital plane 
of the binary and the gas disk can significantly affect 
the dynamics of the planetesimals. Small inclinations 
?r < 10°) can fa vor planetesimal accretion somewhat 
Xie Sz Zhoull2009l ). On the other hand, large inclinations 
(30° < ?b < 50°) can significantly perturb the planetes- 
imal disk, causing planetesimals to "jump" inwards and 
pile up into a smaller inner disk, whe re encounter ve loc- 
ities are more favorable to accretion (|Xie et al.ll2011l ). 

Most of the works in the literature have focused 
on observed or plausible circumstellar configurations 
(e.g. a planet orbiting one of the two stellar com- 
ponents), in light of the lack of direct evidence of 
the existence of circumbinary planets orbiting main- 
sequence stars, outside the realm of science fiction. 
Therefore, only a handful of articles have considered 
planet formation in circumb i nary orbits ( "P-type" ; e.g. , 
[ Moriwaki fc Nakagawall2004t lOuintana fc Lissauerll2006l: 
Scholl et al.ll2007t iMarzari et al.ll2008t iPierens k, Nelson! 
2008ajbl). and they lacked a reference observed configu- 
ration. 

Kepler 16-b (Doyle e t al.l 120111 ) is the first circumbi- 
nary planet that has been detected with Kepler. The 
presence of a third object was first hinted through devi- 
ations of the timing of the stellar eclipses from a linear 
ephemeris. The definitive characterization as a plane- 
tary object came from transits on both star A (tertiary 
eclipse) and star B (quaternary eclipse). The planet 
was determined to be a Saturn-mass planet (A'lp ~ 
0.33A4 j U p) on a nearly circular 228-day orbit; long-term 
integrations have shown the planet to be stable, with an 
eccentricity oscillating between and « 0.08. The binary 
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stellar system is composed of two main-sequence stars in 
an eccentric 41-day orbit, with a mass of 0.69 and 0.2 
M.q (mass ratio fx ~ 0.2), respectively. The close copla- 
narity of the binary and planetary orbital planes suggests 
that the three bodies were formed in a common disk. 
This was bolstered by the me asurement o f the R ossiter- 
McLaughlin doppler shift by iWinn "etaLl (|20ll . which 
indicated t hat the spin of the primary is aligned as well. 

Recently, iWelsh et all (|2012[) reported the discovery of 
two additional circumbinary gas giants (Kepler-34 b and 
Kepler-35 b). The relative abundance of these systems 
among the more than 2,000 ecli psing binaries monitored 
by Kepler (H awson et al.ll20li[) implies a lower limit of 
w 1% in the frequency of circumbinary planets with com- 
parable transit probabilities. Interestingly, all three plan- 
ets lie just outside the stability boundary for test parti- 
cles. Their pericenter distance is, respectively, only « 6% 
(Kepler-34 b), 9% (Kepler-16 b) and 20% (Kepler 35-b) 
larger than the cri tical semi-major axis, as es timated by 
the empirical fit in Holman & Wicgert (1999). This rep- 
resents an important constraint for the formation of the 
planetary core. Indeed, a natural scenario would entail 
the planetary core migrating inwards until near the edge 
of the disk cavity (which will be comparable in extent 
to the stability boundary for test particles), where the 
steep gradient of the disk sur f ace density can hal t migra - 
tion (|Pierens fc Nelsonll2007l ). iPierens fc Nelson! (|2008bn 
simulated the evolution of a 20 A4q core, initially placed 
at the edge of the cavity and free to accrete gas to be- 
come a Saturn-mass planet. They found that once the 
planet depletes the gas in the coorbital region, it will 
resume a slow inward migration, until its eccentricity is 
excited and a phase of runaway outward migration is 
experienced. This runaway migration appeared to stop 
once the planet crossed the 5:1 resonance with the binary, 
at which point slow migration is resumed. The ultimate 
fate of the planet in these simulations is uncertain, due 
to the long timescales involved. However, it is expected 
that disk dispersal will ultimately strand the planet on a 
circular orbit around the binary. Tantalizingly, Kepler- 
16b lies somewhat close (and outside of) the 5:1 period 
ratio with the binary. 

In this paper, we investigate the conditions for the for- 
mation of planetary cores in circumbinary orbits around 
the Kepler-16 binary system, using a simplified numerical 
model. We consider the evolution of a disk of /cm-sized 
planetesimals and determine the impact velocities among 
planetesimals over 10 5 years, the typical timescale fo r 
runaway and oligarchic accretion ([Kokubo &; Idall2000h . 
These preliminary iV-body simulations will be used to 
assess the viability of core accretion as a function of the 
bary centric semi-major axis. 

The plan of the paper is as follows. In JJ21 we briefly dis- 
cuss our numerical model and limitations of our current 
approach. In Sj3]we discuss the results of our simulations 
in the context of planet formation, and conclude in Jll 

2. NUMERICAL SETUP 

To conduct our simulations, we use a new hybrid code, 
Sphiga (described in Meschiari et al., 2012, in prepara- 
tion). Sphiga is an TV-body code that evolves a system 
of non-interacting test particles (e.g. the planetesimals) 
subjected to the sum of gravitational forces of massive 
bodies (e.g. the binary system). In addition, it calculates 



the frictional force acting on the test particles caused by a 
putative protoplanetary disk. By default, this is accom- 
plished by following the complete hydrodynamical evo- 
lution of the disk with the S moothed Particle Hydrody- 
namics scheme (SPH; see, e.g.. !Rosswogl2009t lPricdl2010L 
for recent reviews) in two and three dimensions. The 
same algorithm used to interpolate the hydrodynamical 
quantities can be used to interpolate the local gas density 
and flow and locate possible planetesimal impactors a 
single loop, leading to significant computational savings. 
Modelling the self-consistent perturbations from the bi- 
nary on the disk can alter the planetesi mal evolution 
and p otentially increase impact velocities ()Marzari et al.l 
2008). Indeed, we expect that non-axisymmetric struc- 
ture, such as global spiral patterns, will be imposed by 
the binary, adding a complex time-dependent term. The 
actual impact of the full hydrodynamical evolution is still 
uncertain, however. Even bulk quantities such as the 
disk eccentricity induced by a binary companion appear 
to depend sensitively on the comput ational scheme (e.g. 
the w ave damping prescription in Paardekoopc r et all 
120081 ) and the amount of physics modeled (e.g. the equa- 
tion of state and whether self-gravity was included in 
iMarzari et""aTll2009t IMarzari et al.ll2012t ). 

Nevertheless, significant computational effort is still re- 
quired to follow the evolution of the combined disk, bi- 
nary and planetesimal system (with iV p i + -/V gas > 10 6 
particles) for at least « 10 5 binary revolutions. There- 
fore, for the purpose of this paper, we will use an alter- 
native code path that activates a fixed gas disk. The gas 
disk exerts a frictional acceleration at the location of the 
planetesimal given by 
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In Equation (fTJ), i5v = v p i — v gas is the relative velocity 
of the planetesimal with respect to the Keplerian flow of 
the gas and K is the drag parameter 
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The drag parameter is defined in terms of the planetes- 
imal radius 1Z p i, the planetesimal mass A4 p i (calculated 
assuming p v \ — 3 g/cm 3 ), and the dimensionless coef- 
ficient C (C ~ 0.4 for spherical bodies). We use the 
standard prescription o f a minimum-mass solar nebula 
(MMSN; IHavashi|[l98ll) for the disk parameters. In this 
configuration, our code an d physical setup is fu nctionally 
equivalent to that used bv lScholl et al.l (|2007t ). 

To evaluate the collisional speeds among planetesimals, 
we follow the dynamical evolution of 30,000 test parti- 
cles uniformly distributed with barycentric semi-major 
axes between 0.66 and 6 AU; this range includes the 
current location of the planet (dp « 0.7 AU). The in- 
ner boundary was determined by running a simulation 
with test particles in circular barycentric orbits cover- 
ing semi-major axes in the range (1.2ab;5ab) for 10 4 
years; we found very go od agreement with the fit of 
iHolman fc Wie gert (1999). Particles that travel into the 
inner boundary or become unbound are removed from 
the simulation. 

The system is initially evolved to 10 5 years. After this 
interval, planetesimal-planetesimal close encounters are 
recorded, with the most important parameter being Av, 
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the impac t veloc ity. We follow iFragner et all ([201 1[ ) and 
iThebaultl (|201lD and adopt the prescription for classi- 
fying disruptive impacts for planetesimals presented in 
iStewart fc Leinhardd ((2009). The latter work offers a 
criterion for catastrophic disruption, the main parame- 
ters being the reduced kinetic energy, the masses of the 
impactors and material properties and constants derived 
from fits to numerical and laboratory data. Planetesi- 
mal collisi ons are tracked using the inflated radiu s pre- 
scription (jBrahid [1971 iThebault fc Brahidll999| ). with 
= 5 x 10 -5 AU. The code detects collisions by pop- 
ulating a tree structure at each timestep (as part of the 
SPH algorithm) and walking the tree to locate the near- 
est neighbors to each planetesimal with d < 21Zj„fl (e.g., 
iBarnes fc Huflll986t iHernquist fc Katdll989t ). 

In our simulation, we assign a planetesimal radius for 
each particle, randomly distributed between 1 and 10 
km. We allow for a non-flat primordial distribution in 
planetesimal sizes by assigning a weight /(7^i,7?.2) to 
each impa ct between plan e tesima ls of radius 1Z\ and IZ2 ■ 
Following Pfhebault et all {2008), we use a Maxwellian 
weighting function centered around 5 km with a = 1 
km. A priori, this choice should yield a more accretion- 
friendly environment, since it weighs collisions between 
same-sized planetesimals more than different-sized plan- 
etesimals. 

3. SIMULATIONS 

As expected, the planetesimals are quickly perturbed 
from their initial low-eccentricity configuration by the 
gravitational stirring of the central binary. Their eccen- 
tricities initially oscillate around the forced eccentricity 

e i = \{l-2p)^ (3) 
4 a 

(jMoriwaki fc Nakagawall2004T) . The presence of gas drag 
tends to damp the eccentricity oscillations towards the 
forced eccentricity over longer timescales. Damping and 
periastron phasing will be more effective for smaller 
(since the gas drag coefficient is proportional to 1Z~^ ) and 

close-in bodies (since p gas oc a~ 2 - 75 ). However, the ec- 
centricity spread remains somewhat large at small semi- 
major axes, where the gravitational perturbation of the 
central binary acts to pump eccentricities. At large semi- 
major axes, where the damping timescale is longer, the 
values of eccentricity tend to their counterparts in gas- 
free simulations. 

In the inner parts of the disk, planetesimals will spi- 
ral into the inner boundary due to radial drift. The ra- 
dial drift timescale can be estimated by assuming the 
planetesimal loses angular momentum s lowly due to the 
torqu e from the headwind of the gas ([Weidenschillingl 
1977) . For the drag prescription of Equation^ we find an 
estimate for the infall timescale (in units where GM = 1) 
is given by 

T rd = — ~ 7^ M* — " JTJ (PV) , (4) 

l>rd o Pgas q ^ 

where p v \ is the density of the planetesimal and M.* is 
the total mass of the binary. 

In the case of planet formation around single stars, 
eccentricities are very low and Sv ~ /ipUkep is mainly 
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colorized with respect to their size: light gray (1 < TZ p i < 4 km), 
medium gray (4 < TZ pl < 7 km), black (7 < U pl < 10 km. The 
dashed line shows the forced eccentricity. 

determined by the local scale height ho, with a typical 
drift timescale at 1 AU of 10 6 years for a 5-km plan- 
etesimal. In the circumbinary environment, on the other 
hand, the perturbation from the binary companion acts 
to raise eccentricities throughout the planetesimal disk, 
such that the dominant term contributing to Sv is given 
by the time-varying speed of the planetesimal sampling 
different gas velocities at the apsides. 

Figure [5] shows the distribution of planetesimals after 
t = 10 5 years, binned in semi-major axis. We find that 
inside w 1.5 AU, the planetesimal disk is severely de- 
pleted. Indeed, in our setup, drift timescales are a strong 
function of semi-major axis (oc a~ 5 / 2 ), such that radial 
drift from the outer parts of the disk cannot replenish 
the inner disk effectively. We compared the planetesimal 
distribution of our A^-body run with an analytic model 
based on Equation [3] and @] Assuming Sv w 0.5efWk ep , 
we find good agreement between the two. Finally, the 
second panel of Figure [2] shows that the distribution of 
planetesimal sizes is skewed towards larger planetesimals 
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Semi-major axis (AU) 

Fig. 2. — (Top) Planetesimal number binned in semi-major axis 
after t = 10 s years, normalized by the initial distribution in semi- 
major axis. The TV-body run (black line) and the output from 
an analytic model for the single star case (gray line) and with a 
forced eccentricity term (dashed line) are shown. (Bottom) Rela- 
tive fractions of planetesimals with 1 < 1Z p \ < 4 km (solid line), 
4 < TZ pl < 7 km (dotted line) and 7 < TZ pl < 10 km (dash-dotted 
line) after 10 years. 

at small semi-major axes, since larger planetesimals are 
less affected by the gas drag. This can contribute to mak- 
ing the inner region more accretion-friendly for two rea- 
sons. Firstly, larger planetesimals can withstand larger 
impact velocities. Secondly, the spread in sizes will be 
reduced, which means that the spread in the phasing of 
the planetesimals will also be reduced. 

In the outer parts of the disk, where damping is 
less effective, planetesimals are initially weakly phased 
because the oscillations are coherent and spatially ex- 
tended; therefore, impact velocities tend to be lower. 
However, the frequency of the oscillation around the 
forced eccentricity incre ases with time, ultim ately lead- 
ing to orbital crossing ([Thebault et al.l [20061. The or- 
bital crossing boundary a croS s sweeps outwards in semi- 
major axis, increasing impact velocities. In our simula- 
tion, collisions are recorded for a small time window after 
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Semi-major axis (AU) 

Fig. 3. — Fraction of accreting (dark gray), disturbed (medium 
gray) , uncertain (light gray) and erosive impacts (white) , as a func- 
tion of semi-major axis, after t = 10 5 years. The present-day loca- 
tion of the planet and the fiducial ice line are plotted (dotted and 
dashed line, respectively). 

t = 10 5 years. As evidenced in Figure [TJ regions outside 
« 3.5 AU (w 13ob) have not experienced orbital crossing. 
This is expected, sinc e a cross is a weak function of time 
(|Thebault et al.|[2006f ). Over longer timescales, the im- 
pact velocities will increase in the outer regions, as they 
are swept by the orbital crossing boundary. However, we 
expect the core of Kepler-16 b to be formed and accret- 
ing gas before significant gas dispersal has occurred (on 
a timescale of « 10 5 years; e.g. lWolk fc Walter 1996). 

3.1. Implications for planet formation 

Figure [3] shows the fraction of accreting encounters as a 
function of semi-major axis. We found that the following 
qualitative situation holds for different radial locations: 

(a) in the region between the stability boundary and 1 
AU (which includes the present-day location of the 
planet ap ~ 0.7 AU), eccentricities are pumped to 
high values by the central binary and planetesimal 
number density is low due to the fast radial drift. 
The majority of encounters are in the "uncertain" 
regime, with the potential of being accreting depend- 
ing on the prescription for the outcome of disruptive 
collisions. 

(b) for a small range in semi-major axis outside 1 AU, 
the spread in e and w is smaller and planetesimal dis- 
tributions are skewed towards larger planetesimals. 
The majority of encounters are accreting. 

(c) between 1.75 AU and 4 AUs, the magnitude of the 
eccentricity and the differential phasing raises the im- 
pact velocities, such that the majority of the encoun- 
ters are erosive. 

(d) outside 4 AUs, orbital crossing has not been realized 
yet and gas drag is weaker due to the steep radial 
dependence of the gas density; therefore, orbits are 
only weakly phased. The majority of encounters are 
accreting. 

We conclude that planet formation is likely inhibited 
for a large range in semi-major axis (location (c), between 
1.75 and 4 AUs). This range in semi- major axis includes 
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Fig. 4. — Initial location of embryos that grow to final masses 
0.2 < Mp < 0.4.Mj (black points) and Mp _ n3l i > lMj (red dots) 
in the simulations of Mordasini ct al. (2012]), for a range of metal- 
licities and disk masses. The shaded region corresponds to the 
range in semi-major axis where embryo formation is disturbed in 
the Kepler-16 system. 

the nominal location of the ice line for an irradiated disk, 
estimate d from the scal ing: a ice ~ 2.7AU (M/Mq) 2 « 
2.3 AU (llda fc Linl feOtM). assuming M = Ai\ + A^b- 

What is the impact of this "forbidden region" for 
planet formation? It is instructive to refer to the pre- 
dictions of the standard core accretion paradigm for 
single stars; in particular, the outcome of large-scale 
Monte-Carlo planet synthe s is models (e.g. llda fc Linl 
[200l iMordasini et al.ll2009h . iMordasini et all (|2012t ) r~ 
cently conducted a Monte-Carlo planet synthesis simula- 
tion for a variety of disk masses and metallicities, for the 
nominal case of a 1 .A/f© central star. In the core accre- 
tion paradigm, metallicity represents a threshold quan- 
tity for the formation of planetary cores. Accordingly, 
they found that the cores of giant planets (M > Mj) 
tend to preferentially form outside the ice line when the 
metallicity (which acts as a proxy for the solid content of 
the disk) is low. The actual location of the ice line scales 
with the disk mass, which contributes to the spread in 
semi-major axis. 

In Figure (g]), we p lot a different subset of the output 
of the simulations of IMordasini et al~l ([20120, focusing 
on the ensemble of embryos that acquire masses com- 
parable to Kepler-16 b (0.2Mj < M < OAMj). The 
initial location of the embryo is plotted as a function of 
metallicity. For disks of solar or super-solar metallicity, 
such planets are formed throughout the disk, with a sub- 
stantial fraction formed inside 2 AU (about 40%). At 
subsolar metallicities comparable to Kepler-16 ([Fe/H] 
w — 0.3 ± 0.2), however, such cores are only found out- 
side 2 AU, with a minority lying in location (c) (about 
20%). While the synthetic population refers to the nom- 
inal 1 M.q single star case, with one embryo per disk, 
it suggests that in situ planet formation in location (a) 
might be hampered by the low surface density in solids 
at 1 AU. Our simulations place an additional dynamical 
constraint, indicating that less than 20% of encounters 

3 http:/ /www. mpia-hd.mpg.de/homes/mordasini/Site7. html 



within 1 AU are accreting. This, compounded with the 
low planetesimal density in the region (Figure [2]) , makes 
in situ formation of a substantial core difficult. 

Finally, it is also crucial to recognize that non- 
axisymmetric perturbations from the disk might play an 
important role in the dynamics of the inner disk. The ec- 
centric central binary will likely excite spiral structures, 
which might act to pump the eccentricity of the inner 
planc tesimals and a l ter th e phasing of their orbits. In- 
deed, iMarzari et all (|2008h conducted full 2D hydrody- 
namical simulations with a small number of tracer plan- 
etesimals embedded in the disk, and found significant os- 
cillations in the eccentricity and longitude of pericenter 
around the equilibrium value. 

4. DISCUSSION 

Planet formation in presence of close binaries presents 
a number of challenges to the traditional core accretion 
paradigm. Historically, most of the theoretical effort 
has been expended to study pathways to planet forma- 
tion in S-type orbits for planets that had been observed 
through RV surveys, or targets with observationally de- 
sirable properties (e.g., a Centauri). With the launch of 
Kepler, however, we expect that the sample of planets in 
P-type orbits around eclipsing binaries will rapidly out- 
number the handful of planets in circumstellar configura- 
tions detected with RV surveys. Indeed, a sample of 750 
Kepler targets are eclipsing binaries for which eclipses 
of both stars are observed, and a subset o f 18% exhib- 
ited d eviations in the timing of the eclipses ([Welsh et ID 
Since the definitive determination of the plane- 
tary nature of a putative KOI relies on the detection of 
tertiary and quaternary eclipses, we expect that as the 
baseline of the observation increases, more KOIs will be 
confirmed as genuine circumbinary objects. 

In this paper, we have conducted a preliminary simula- 
tion of the feasibility of circumbinary planet formation in 
the Kepler-16 system. In accord ance to an earlier study 
conducted bv lScholl et a l. (2007) for a different set of bi- 
nary parameters, we have found that, for generous initial 
conditions that favor planetesimal accretion, planet for- 
mation appears to be feasible far enough from the cen- 
tral binary. However, we have identified a substantial 
radial span between 1.75 and 4 AU where planet forma- 
tion is strongly inhibited. Within the planet accretion 
framework, the most likely sequence of event is the for- 
mation of a core outside the forbidden region, followed 
by inwards migration driven by tidal interact ion with the 
protoplanetary disk ()Pierens fc Nelsonll2007l) . Although 
we measured impact velocities potentially favorable to 
accretion close to the present-day location of the planet, 
in situ formation of Kepler-16 b is less likely due to over- 
all high encounter speeds, low planetesimal density, low 
metallicity of the star, and non-axisymmetric perturba- 
tions from the disk (not modeled in this simulation). 

We remark that the simulations presented in this paper 
only demonstrate that, choosing the most favorable con- 
ditions for planetesimal accretion and an assumed initial 
planetesimal size of 1-10 km, the formation of an embryo 
outside 4 AU is plausible, with traditional migration pro- 
cesses subsequently moving the planet to its current loca- 
tion. Our approach has several limitations introduced for 
the sake of simplicity and computational speed; chiefly, 
we disregarded the evolution of the protoplanetary disk 
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and the collisional outcome of planetesimal impacts. For 
the former, we plan to follow approximately the hydro- 
dynamical response of the disk with the SPH algorithm 
included in the Sphiga code in a follow-up paper. For 
the latter, a time-dependent distribution of planetesi- 
mal sizes would more accurately model the extent of the 
accretion-friendly regions, which depend sensitively on 
the planetesimal parameters. The numerical procedure 



of iPaardekooper fc Leinhardt] (|2010ft represents a possi- 
ble approach to following the collisional evolution of the 
planetesimal size distribution. 

We acknowledge support from the NASA Grant 
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